Abstract
Background In the report naomi-simple_fit with parameter tmbstan = TRUE, we used the NUTS algorithm to perform MCMC inference for the simplified Naomi model.
Task Here we assess whether or not the results of the MCMC are suitable using a range of diagnostic tools.
We start by obtaining results from the latest version of naomi-simple_fit with tmbstan = TRUE.
out <- readRDS("depends/out.rds")
mcmc <- out$mcmc$stanfit
This MCMC took 3.28 days to run
cbpalette <- c("#56B4E9", "#009E73", "#E69F00", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999")
bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())
We are looking for values of \(\hat R\) less than 1.1 here.
rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat_data(rhats) %>%
ggplot(aes(x = value, y = parameter, col = description)) +
geom_segment(aes(yend = parameter, xend = ifelse(min(value) < 1, 1, -Inf)), na.rm = TRUE, alpha = 0.7) +
scale_color_manual(values = "#E69F00") +
geom_vline(xintercept = 1.05, linetype = "dashed", col = "grey40") +
labs(x = "Potential scale reduction factor", y = "NUTS parameter", col = "") +
theme_minimal() +
theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank()
)
ggsave("rhat.png", h = 3, w = 6.25)
(big_rhats <- rhats[rhats > 1.1])
## named numeric(0)
length(big_rhats) / length(rhats)
## [1] 0
Reasonable to be worried about values less than 0.1 here.
ratios <- bayesplot::neff_ratio(mcmc)
breaks <- c(0, 0.1, 0.25, 0.5, 0.75, 1)
bayesplot::mcmc_neff_data(ratios) %>%
ggplot(mapping = aes(x = value, y = parameter, color = description)) +
geom_segment(aes(yend = parameter, xend = -Inf), na.rm = TRUE, alpha = 0.7) +
scale_color_manual(values = c("#56B4E9", "#009E73", "#E69F00")) +
geom_vline(xintercept = 0.1, linetype = "dashed", col = "grey40") +
geom_vline(xintercept = 0.5, linetype = "dashed", col = "grey40") +
geom_vline(xintercept = 1, linetype = "dashed", col = "grey40") +
labs(x = "ESS ratio", y = "NUTS parameter", col = "") +
theme_minimal() +
theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank()
)
ggsave("ratio.png", h = 3, w = 6.25)
What are the total effective sample sizes?
#' I think that this $summary should be all of the chains grouped together
mcmc_summary <- summary(mcmc)$summary
data.frame(mcmc_summary) %>%
tibble::rownames_to_column("param") %>%
ggplot(aes(x = n_eff)) +
geom_histogram(alpha = 0.8) +
labs(x = "ESS", y = "Count")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("ess.png", h = 3, w = 6.25)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
How much autocorrelation is there in the chains?
bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model
Variation between units can be explained by high correlation and high standard deviation, or low correlation and low standard deviation. Hence there is an unidentifiabiility here that leads to correlated posteriors.
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_alpha_as", "logit_phi_alpha_as"), diag_fun = "hist", off_diag_fun = "hex")
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_rho_as", "logit_phi_rho_as"), diag_fun = "hist", off_diag_fun = "hex")
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_alpha_a", "logit_phi_alpha_a"), diag_fun = "hist", off_diag_fun = "hex")
(plot <- bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_rho_a", "logit_phi_rho_a"), diag_fun = "hist", off_diag_fun = "hex"))
ggsave("rho_a.png", plot, h = 4, w = 6.25)
There is a prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable.
Let’s have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter.
area_merged <- sf::read_sf(system.file("extdata/demo_areas.geojson", package = "naomi"))
nb <- area_merged %>%
filter(area_level == max(area_level)) %>%
bsae::sf_to_nb()
neighbours_pairs_plot <- function(par, i) {
neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}
# area_merged %>%
# filter(area_level == max(area_level)) %>%
# print(n = Inf)
Here are Nkhata Bay and neighbours:
neighbours_pairs_plot("log_or_gamma", 5)
And here are Blantyre and neighbours:
neighbours_pairs_plot("log_or_gamma", 26)
np <- bayesplot::nuts_params(mcmc)
saveRDS(np, "nuts-params.rds")
Are there any divergent transitions?
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_divergent = sum(Value))
bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))
We can also use energy plots (Betancourt 2017): ideally these two histograms would be the same When the histograms are quite different, it may suggest the chains are not fully exploring the tails of the target distribution.
bayesplot::mcmc_nuts_energy(np)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 13.3.1
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] multi.utils_0.1.0 sf_1.0-9 bayesplot_1.9.0 patchwork_1.1.2 inf.utils_0.1.0 aghq_0.4.1 tmbstan_1.0.4 rstan_2.21.5 StanHeaders_2.21.0-7 TMB_1.9.2
## [11] Matrix_1.5-4.1 forcats_0.5.2 stringr_1.5.0 dplyr_1.0.10 purrr_1.0.1 readr_2.1.3 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.0 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 deldir_1.0-6 ellipsis_0.3.2 class_7.3-20 ggridges_0.5.3 rprojroot_2.0.3 fs_1.6.1 rstudioapi_0.14 proxy_0.4-27 hexbin_1.28.2
## [11] farver_2.1.1 fansi_1.0.4 mvtnorm_1.1-3 lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18 splines_4.2.0 cachem_1.0.6 knitr_1.41 polynom_1.4-1
## [21] jsonlite_1.8.4 broom_0.8.0 naomi_2.8.5 dbplyr_2.1.1 compiler_4.2.0 httr_1.4.4 backports_1.4.1 assertthat_0.2.1 fastmap_1.1.0 cli_3.6.1
## [31] s2_1.1.1 htmltools_0.5.3 prettyunits_1.1.1 tools_4.2.0 gtable_0.3.1 glue_1.6.2 reshape2_1.4.4 posterior_1.2.2 wk_0.7.0 V8_4.2.2
## [41] Rcpp_1.0.10 jquerylib_0.1.4 cellranger_1.1.0 vctrs_0.6.1 spdep_1.2-7 mvQuad_1.0-6 tensorA_0.36.2 xfun_0.37 ps_1.7.3 rvest_1.0.2
## [51] bsae_0.2.7 lifecycle_1.0.3 statmod_1.4.36 orderly_1.4.3 ids_1.0.1 MASS_7.3-56 scales_1.2.1 ragg_1.2.2 hms_1.1.2 parallel_4.2.0
## [61] inline_0.3.19 yaml_2.3.7 curl_5.0.0 gridExtra_2.3 sass_0.4.4 loo_2.5.1 stringi_1.7.8 highr_0.9 checkmate_2.1.0 e1071_1.7-12
## [71] boot_1.3-28 pkgbuild_1.3.1 spData_2.2.1 rlang_1.1.0 pkgconfig_2.0.3 systemfonts_1.0.4 matrixStats_0.62.0 distributional_0.3.0 evaluate_0.20 lattice_0.20-45
## [81] labeling_0.4.2 cowplot_1.1.1 processx_3.8.0 tidyselect_1.2.0 traduire_0.0.6 bookdown_0.26 plyr_1.8.8 magrittr_2.0.3 R6_2.5.1 generics_0.1.3
## [91] DBI_1.1.3 pillar_1.9.0 haven_2.5.0 withr_2.5.0 units_0.8-0 sp_1.5-1 abind_1.4-5 modelr_0.1.8 crayon_1.5.2 uuid_1.1-0
## [101] KernSmooth_2.23-20 utf8_1.2.3 tzdb_0.3.0 rmarkdown_2.18 grid_4.2.0 readxl_1.4.0 isoband_0.2.6 data.table_1.14.6 callr_3.7.3 reprex_2.0.1
## [111] digest_0.6.31 classInt_0.4-8 numDeriv_2016.8-1.1 textshaping_0.3.6 openssl_2.0.5 RcppParallel_5.1.5 stats4_4.2.0 munsell_0.5.0 bslib_0.4.1 askpass_1.1